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Nickel  zirconia,  Ni-YSZ,  porous  cermets  used  in  solid  oxide  fuel  cells  can  be  considered  as  composite 
materials  consisting  of  three  phases,  namely  nickel  (Ni),  ytrria-stabilized  zirconia  (YSZ),  and  voids.  Based 
on  Ni-YSZ’s  microstructural  features,  such  as  the  volume  fractions,  average  particle  sizes  and  their  cor¬ 
relations  of  a  physical  sample,  a  three-dimensional  stochastic  reconstruction  method  is  used  to  recreate 
the  three-phase  composite.  The  digitally  reconstructed  microstructure  is  then  transferred  into  finite  ele¬ 
ment  models  to  obtain  the  material’s  corresponding  effective  elastic  modulus  and  effective  coefficient  of 
thermal  expansion  at  various  temperatures.  Predictions  from  such  a  numerical  methodology  agree  well 
with  experimental  results. 

©  2008  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Planar  solid  oxide  fuel  cells  (SOFCs)  are  high-density  power 
sources  which  promise  relatively  small  environmental  impact.  Yet 
their  widespread  use  is  hindered  by  high  costs  and  poor  long-term 
thermomechanical  reliability.  Part  of  this  reliability  issue  extends 
from  thermal  mismatch  between  the  fuel  cell  layers  and  strict 
mechanical  constraints  from  the  SOFC  configuration.  For  example, 
SOFCs  are  constructed  of  repeating  layers  of  porous  composites  and 
solid  ceramics  across  which  the  electrochemical  reactions  will  take 
place.  The  composites  are  constructed  of  interpenetrating  metals 
and  ceramics,  called  cermets,  that  provide  structural  stability,  effi¬ 
cient  flow  paths  for  fuel  and  air,  and  conduction  pathways  for  ion 
transfer.  Each  of  these  functions  is  influenced  by  the  respective 
arrangement  and  amount  of  metal  and  ceramic.  Therefore  there 
is  a  need  to  understand  how  the  cermet’s  microstructure  influ¬ 
ences  the  thermomechanical  properties  of  interpenetrating  metal 
ceramic  composites. 

This  work  numerically  creates  multiple  realizations  (computer 
generated  digital  models  of  the  material)  of  these  composites  to 
serve  as  a  platform  for  numerical  analyses.  The  realizations  can 
then  be  used  to  determine  material  properties  and  their  links  to 
the  microstructure. 
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The  anode  material,  nickel-zirconia  (Ni-YSZ),  is  a  classic  exam¬ 
ple  of  a  three-phase  porous  material  with  a  random  heterogenous 
microstructure.  Made  by  tape  casting  NiO-YSZ  slurry,  the  final 
microstructure  is  highly  influenced  by  starting  particle  sizes,  con¬ 
stituent  compositions,  and  the  multi-step  manufacturing  process. 
The  final  microstructure  is  also  influenced  by  the  oxidation  of 
NiO  during  the  sintering  process  and  fuel  cell  operation.  Although 
Ni-YSZ’s  primary  use  is  as  a  vehicle  for  ion  transfer,  the  material 
must  also  meet  several  other  requirements.  It  must  provide  struc¬ 
tural  support  for  the  other  layers  in  the  cell;  it  must  allow  sufficient 
flow  paths  for  fuels,  and  finally  it  must  be  able  to  handle  tempera¬ 
ture  cycles  from  room  temperature  up  to  several  hundred  degrees 
Celsius.  These  requirements  placed  on  the  anode  require  a  wide 
range  of  different  material  behaviors,  and  the  enhancement  of  any 
of  these  can  negatively  impact  another  property. 

In  this  paper  Ni-YSZ  is  numerically  generated  as  a  three-phase 
composite  made  of  overlapping  spheres  of  nickel  and  YSZ  with  a 
porous  matrix.  The  digital  reconstruction  method  recreated  a  three- 
dimensional  pixilated  image,  or  realization,  of  the  material.  Since 
the  reconstruction  is  3D,  the  term  voxel  will  be  used  in  place  of 
pixel.  By  generating  multiple  realizations,  the  statistical  variation 
of  key  material  properties,  such  as  modulus  and  the  coefficient  of 
thermal  expansion  (CTE),  can  be  determined. 

Reconstruction  and  analysis  of  heterogenous  structures  takes 
place  down  several  avenues.  Fiber  and  particulate  composites  can 
be  reconstructed  through  tessellation  procedures  and  then  ana¬ 
lyzed  using  self-consistent  methods.  Pyrz  studied  fiber  composites 
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using  Dirichelet  tessellations  for  simulated  hard  core  models  and 
microscopy  images,  respectively  [1  ].  Further  research  by  Bochenek 
and  Pyrz  studied  unidirectional  fiber  reinforced  composites  and 
particulate  composites  with  Voronoi  tessellations  [2].  Stress  inter¬ 
actions  around  the  particles  were  calculated  using  superposition 
and  the  effective  medium  approach. 

Ghosh  and  co-workers  have  directly  incorporated  representa¬ 
tive  volume  elements  (RVEs)  generated  from  Voronoi  cell  (VC) 
tessellations  into  a  multi-scale  finite  element  method  (FEM)  [3-6]. 
The  microstructural  model  is  termed  VC-FEM,  and  it  models  par¬ 
ticulates  in  a  matrix  using  Voronoi  tessellations,  where  each 
tessellation  is  treated  as  an  element  with  a  particle  at  its  center  [3]. 
The  same  method  was  used  for  a  multi-scale  damage  analysis  in 
porous  materials  [5].  Additionally,  three-dimensional  models  were 
created  through  stereological  methods  for  particle  reinforced  metal 
matrix  composites  [4].  Next  the  concept  of  statistically  equivalent 
representative  volume  elements  (SERVEs)  was  used  with  VC-FEM 
to  study  fiber  and  particulate  composites  that  varied  randomly 
within  the  microstructure  [7,8].  Methods  developed  by  Gokhale 
and  co-workers  used  extremely  detailed  images  of  the  particulates 
in  heterogenous  microstructures  to  recreate  the  material  through 
a  controlled  re-distribution  methodology  [9-12].  Another  work 
examined  the  most  efficient  RVE  size  for  use  in  macro-micro  analy¬ 
sis  for  fiber  composites  [13].  Other  methods  perform  finite  element 
analyses  on  the  actual  microstructures  read  in  via  open  source  soft¬ 
ware  OOF;  then  the  FE  model  can  be  used  for  different  studies,  such 
as  Cannillo’s  and  Carter’s  stochastic  damage  analysis  on  a  polycrys¬ 
talline  microstructure  [14]. 

A  modification  of  the  simulated  annealing  method  is  used  here 
for  its  flexibility  and  efficiency  in  creating  multiple  realizations 
for  numerical  analysis.  Rintoul  and  Torquato  used  the  simulated 
annealing  method  to  reconstruct  a  distribution  of  spheres  using 
the  radial  distribution  function  [15].  Next  the  simulated  anneal¬ 
ing  method  was  used  to  recreate  digitized  media  by  Yeong  and 
Torquato  [16].  Yeong  and  Torquato  minimized  processing  time  by 
calculating  the  2-point  correlation  function  and  the  lineal  chord 
function  in  orthogonal  directions  and  then  updating  the  functions 
only  along  rows  and  columns  with  one  to  one  pixel  exchanges.  Man- 
wart  and  Hilfer  noted  that  the  time  saving  device  of  calculating 
correlation  functions  in  only  orthogonal  directions  will  only  intro¬ 
duce  anisotropy  effects  for  microstructures  displaying  significant 
short-range  order  [17].  The  impact  of  short-range  order  and  the 
behavior  of  the  simulated  annealing  procedure  is  further  analyzed 
by  Cule  and  Torquato  [18].  Microstructural  information  from  a  2D 
slice  can  also  be  used  to  construct  a  three-dimensional  image  as 
shown  in  Part  II  of  Yeong’s  and  Torquato’s  work  on  reconstruct¬ 
ing  random  media  [19].  Rozman  and  Utz  used  several  techniques 
to  improve  the  efficiency  of  the  Monte  Carlo  reconstruction;  these 
include  using  the  Great  Deluge  Algorithm  plus  an  additional  crite¬ 
rion  for  “uphill”  moves,  limiting  pixel  changes  to  the  interface,  and 
calculating  perturbations  of  the  probability  functions  [20]. 

The  material  properties  for  random  media  can  be  determined 
through  analytical,  numerical,  or  experimental  means.  However, 
each  method  has  limitations.  For  random  co-continuous  media 
analytical  approaches  are  limited  by  microstructural  information, 
and  this  is  especially  true  for  porous  media.  Unit  cell  models 
will  be  limited  to  very  simple  composites.  Variational  methods  to 
determine  bounds,  such  as  Hashin-Sthrikman,  will  diverge  away 
from  the  upper  bound  due  to  porosity.  Experimental  methods, 
while  vital,  are  also  severely  limited  by  expense  and  efficiency. 
While  software  such  as  OOF  can  input  one  image  of  a  material, 
it  becomes  significantly  more  difficult  to  study  the  variation  of 
relevant  material  properties.  Additionally,  two-dimensional  recon¬ 
struction  techniques  such  as  tessellations  can  accurately  capture 
the  behavior  of  a  material,  but  cannot  be  easily  adapted  to  the 


behavior  of  a  three-dimensional  random  interpenetrating  compos¬ 
ite.  To  that  end  numerical  reconstruction  methods  in  connection 
with  finite  element  analysis  can  become  a  fast  and  flexible  way  to 
determine  realistic  material  properties  with  specific  microstruc¬ 
tures. 

The  simulated  annealing  method  will  be  used  to  create  voxel 
reconstructions  for  numerical  analysis.  The  effective  modulus  and 
coefficient  of  thermal  expansion  will  be  found  using  multiple  FE 
analyses.  These  analyses  will  also  investigate  size  effects  of  the 
voxel  reconstruction  and  the  statistical  variation  of  the  results. 

The  following  work  will  be  arranged  as  follows.  First,  Sec¬ 
tion  2  will  introduce  the  methodology  behind  random  media  and 
probability  functions.  Section  2.3  will  specifically  cover  the  recon¬ 
struction  method  used,  and  Section  2.4  will  introduce  the  FE  model. 
Section  3  will  go  over  the  results  of  the  reconstruction,  model  con¬ 
vergence,  and  the  effective  properties.  The  last  two  sections  will 
discuss  these  results  in  detail. 


2.  Methodology 

2.1.  Random  media 


For  any  media  of  volume  Vz  the  microstructure  can  be  fully  char¬ 
acterized  by  an  indicator  function 


if  x  e  Vt 
otherwise  ’ 


(2.1) 


where  i  is  the  phase  number,  and  x  is  the  vector  location  within 
the  volume  [21].  The  microstructure  is  assumed  to  be  static  and 
therefore  is  not  a  function  of  time. 

The  indicator  function  describes  every  possible  point  with  a 
material,  such  that  for  any  number,  k,  of  phases  the  following  equal¬ 
ity  holds, 


k 

£/«(*)  =  i. 

i=l 


(2.2) 


With  the  indicator  function  random  media  can  be  described  by 
determining  the  probability  of  a  desired  event  or  occurrence.  For 
example,  an  event  of  interest  could  be  when  multiple  points  lie 
within  the  same  phase.  Such  an  event  is  an  example  of  the  n- point 
probability  function  as  illustrated  in  (2.3). 

S®(x,,  x2,  . . . ,  x„)  =  P(J®(x i)  =  l,/®(x2)  =  l,  ...,  J®(x„)  =  l), 

(2.3) 

where  P  indicates  the  probability  that  a  given  location  lies  within 
phase  i.  As  the  order  n  increases  more  microstructural  details  are 
captured. 

If  the  probability  distributions  of  a  material  are  invariant  with 
respect  to  location,  the  material  is  statistically  homogenous  and  the 
material  is  called  ergodic.  Additionally,  if  the  random  media  does 
not  depend  on  the  orientation  of  the  vector  positions,  but  only  on 
the  magnitude  of  the  distance  between  the  points,  it  can  be  con¬ 
sidered  isotropic.  In  this  case,  the  n-point  functions  now  become 
functions  of  the  distance  between  the  points  such  that  r  =  \Xj  -  x,-|. 

Thus  for  homogenous  media  the  1 -point  function  will  reduce  to 
the  volume  fraction  [cpO  of  the  material, 

S®  =  P(J®(x,))  =  (Pi,  (2.4) 

and  the  2-point  functions  become  functions  of  distance  r, 

S®(r)  =  S®(x1;  x2)  =  P(/®(X!)  =  1,  J®(x2)  =  1).  (2.5) 
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Bounds  exist  for  the  2-point  function  in  homogenous  media  as 
the  radius  reduces  to  zero  or  extends  to  infinity.  These  are 

limS®(r)  =  ^  (2.6) 

r-^0 

and 

limS'i)(r)  =  te)1 2-  (2.7) 

r^oo 

In  Eq.  (2.6)  the  2-point  functions  reduce  to  the  1 -point  function 
as  r  decreases  and  the  two  points  converge  to  each  other.  In  Eq.  (2.7), 
as  the  distances  between  the  points  increase  they  are  no  longer 
spatially  correlated  and  the  2-point  approach  becomes  equivalent 
to  calculating  the  1 -point  function  at  two  separate  points. 

The  behavior  of  the  2-point  probability  function  as  r  approaches 
zero  is  the  first  hint  of  the  relationship  between  the  probability 
functions  for  multiple  phases.  Torquato  and  Stell  showed  that  any  n- 
point  probability  function  can  be  written  as  a  function  of  the  other 
phases  [22].  In  other  words,  for  a  two-phase  material  the  descrip¬ 
tion  of  one  phase  will  guarantee  its  complement  to  the  second 
phase.  However,  as  the  number  of  phases  increases  the  relationship 
must  also  be  quantified  between  the  phases,  i.e.  for  three  phases 
there  are  actually  multiple  phases,  namely  phase  1,  phase  2,  phase 
3,  and  the  combinations  of  any  two  phases. 

2.2.  Microstructure 

For  this  work  the  nickel  and  zirconia  phases  are  described  as 
overlapping  spheres  as  defined  by  Weissberg  [21,23].  The  2-point 
probability  function  does  not  exhibit  any  short-range  order  and  is 
defined  by  one  characteristic  length,  the  particle  diameter  (d).  At 
this  point  the  function  has  reached  its  long-range  order  shown  in 

(2.7) .  While  the  function  controls  particle  size  it  does  not  corre¬ 
late  the  particles  together,  which  allows  interpenetration  to  occur, 
and  the  most  significant  factor  on  clustering  will  be  the  volume 
fraction.  The  analytical  expression  for  overlapping  spheres  uses  the 
Heaviside  function  (0)  as  shown  in  (2.8). 

S®  =  1  -  2(1  -  (Pi)  +  (1  -  <pi)(4s(r)/5rd2), 

where g(r)  =  (:r  -  ®(d  -  r)(arccos(r/d)  -  (r/d)^/l  -  r/d)) 

(2.8) 

A  third  function  is  now  required  to  create  each  realization.  For 
the  three  phases,  the  porous  phase  was  the  matrix  of  the  composite ; 
however,  since  no  analytical  function  exists  for  this  value  it  was 
approximated  from  one  run  of  the  material.  The  estimated  2-point 
function  for  the  pore  phases  will  satisfy  the  bounds  in  (2.6)  and 

(2.7) ,  but  can  exhibit  short-range  order. 

2.3.  Reconstruction 

The  realizations  were  generated  using  the  digitized  simulated 
method  introduced  by  Lee  and  Torquato  with  modifications  from 
Rozman  and  Utz  [3,20].  The  algorithms  were  implemented  in  the 
C++  language  and  the  GNU  GCC  compiler  [24].  The  reconstruction 
procedure  modifies  the  indicator  function  in  Eq.  (2.2)  until  the  sam¬ 
ple  matches  the  desired  probability  functions.  Then  the  indicator 
function  is  used  to  create  a  voxel  representation  of  the  material  for 
further  numerical  use.  Each  voxel  represents  a  different  phase.  The 
multiple  step  process  is  described  below: 

1.  Randomly  seed  sample  with  volume  fractions  of  each  materials. 

2.  Calculate  desired  probability  functions  in  the  current  realization. 


Fig.  1.  2-Point  probability  functions  for  reconstructed  Ni-YSZ  where  nickel  and  YSZ 
are  overlapping  spheres  with  a  5  p,m  diameter. 

3.  Exchange  random  voxels.  The  exchanging  of  voxels  was  con¬ 
strained  such  that  the  volume  fraction  of  all  three  phases  was 
maintained. 

4.  Update  probability  functions  of  volume  and  calculate  energy. 

5.  Accept  or  reject  change. 

6.  Repeat  until  an  accepted  criterion  is  reached. 

The  efficiency  of  the  process  was  increased  by  using  voxel  selec¬ 
tion  at  the  interface,  by  sampling  in  orthogonal  directions,  by  using 
the  great  deluge  algorithm,  and  by  using  perturbations  to  calculate 
the  correlation  functions.  Boundaries  are  periodic. 

Part  of  the  flexibility  of  the  DIB  reconstruction  is  the  flexibility 
inherent  in  the  use  of  energy  as  an  acceptance  criterion  [16].  Any 
desired  correlation  function  can  be  used  and  these  individual  func¬ 
tions  can  be  weighted  as  desired.  The  energy  function  minimizes 
the  least  square  difference  between  a  given  set  of  desired  proba¬ 
bility  functions  and  the  functions  existing  at  the  current  time  step 
of  the  reconstruction.  Use  of  the  least  squares  difference  allows 
flexibility  in  the  reconstruction  by  allowing  different  probability 
functions  to  be  used  at  different  weights. 

2.4.  FE  model 

The  commercial  software  ABAQUS  was  used  to  perform  the  finite 
element  analysis.  For  each  digital  image  reconstructed,  the  voxels 
were  transferred  to  eight  node  brick  elements  and  input  into  a  FE 
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Table  1 

Constituent  material  properties  used  in  FE  analysis  at  25  °C 


model.  Perfect  bonding  was  assumed  between  the  elements,  and 
each  element  was  assigned  the  material  properties  corresponding 
to  the  digital  reconstruction.  The  Young’s  modulus  and  coefficient 
of  thermal  expansion  of  the  volume  were  calculated  in  each  coor¬ 
dinate  direction.  The  Young’s  modulus  was  calculated  at  room 
temperature  using  a  standard  structural  analysis,  while  the  CTE  was 
calculated  for  a  temperature  range  between  20°  and  1000°. 

The  constituent  material  properties  were  temperature  depen¬ 
dant  and  were  taken  from  the  literature  referenced  in  Table  1. 
Experimental  studies  of  the  electrolyte  material,  non-porous 
8  mol%  ytrria  (YSZ)  material,  were  used  as  the  YSZ  properties  in  the 
Ni-YSZ  composite.  The  behavior  of  nickel  around  the  Curie  point 
is  complicated  by  the  sudden  jump  in  thermal  expansion  from  the 
ferromagnetic  transition;  therefore  a  temperature  dependent  equa¬ 
tion  of  CTE  was  used  from  the  work  of  Faisst  [25].  This  expression 
was  used  to  provide  the  CTE  from  0  to  1000  °C.  The  CTE  values  for 
nickel  and  YSZ  are  included  in  Fig.  9. 

3.  Results 

3.1.  Realization  reconstruction 

Ni-YSZ  is  usually  made  from  a  precursor  of  NiO-YSZ,  and  a 
commonly  used  composition  in  fuel  cells  is  75%  NiO/25mol%  YSZ 
[26].  When  the  anode  material  is  fully  reduced  to  Ni-YSZ,  stud¬ 
ies  have  found  the  volume  porosity  to  be  40%  and  the  nickel  and 
YSZ  ratios  can  be  calculated  to  be  25%  and  35%,  respectively  [26]. 
For  both  nickel  and  YSZ,  the  sphere  diameter  was  set  to  5  pm 
and  the  initial  total  length  of  the  sample  was  50  pm.  One  realiza¬ 
tion  was  generated  with  nickel  and  YSZ  modeled  as  overlapping 
spheres.  The  pore  probability  function  from  this  realization  was 
smoothed  for  long-range  order  and  then  used  for  all  further  real¬ 
izations.  Fig.  1  shows  the  2-point  probability  function  for  all  three 
phases  for  a  realization  generated  from  all  three  probability  func¬ 
tions.  The  pore  function  shows  short-range  order  that  vanishes 
before  8  pm.  Fig.  2  compares  the  2-point  probability  function  of 
porosity  to  the  combined  nickel  and  zirconia  phases.  It  can  be  seen 
that  the  length  and  magnitude  of  the  short-range  order  are  com¬ 
parable.  The  2-point  pore  function  plotted  in  Fig.  1  is  used  for  all 
future  realizations  and  now  each  phase  is  governed  by  a  probability 
function. 

An  additional  probability  function  was  calculated  from  the 
finished  realization.  The  lineal  chord  function  measures  the  prob¬ 
ability  that  a  chord  connecting  two  random  points  in  a  phase  i  will 
lie  in  the  same  phase  i.  The  lineal  chord  functions  for  each  phase 
are  plotted  in  Fig.  3.  Examination  of  the  plot  shows  that  no  signifi- 
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Fig.  3.  Lineal  chord  functions  for  each  phase  in  the  realization. 


cant  chord  lengths  exist  past  15  pm,  approximately  three  times  the 
particle  diameters  of  nickel  and  YSZ. 

3.2.  Model  convergence 

After  converting  the  realization  to  the  FE  model,  convergence 
was  tested  by  two  different  methods:  discretization  error  and 
representative  volume  element  size.  Mesh  density  was  tested  by 
subdividing  each  voxel  into  four  subsections  and  rerunning  the 
analysis,  and  no  significant  difference  was  found  to  occur.  Exam¬ 
ination  of  discretization  error  determines  that  the  voxel  size  is 
small  enough  to  accurately  capture  material  behavior,  while  the 
RVE  size  ensures  that  a  large  enough  sample  of  the  material  is 
analyzed  to  capture  the  material  behavior.  Separate  tests  were 
performed  for  both  the  modulus  and  the  CTE,  since  RVE  size 
and  discretization  error  could  vary  depending  on  the  analysis 
type. 

3.2.1.  Discretization 

Table  2  presents  the  calculated  modulus  for  five  different  real¬ 
izations  in  three  directions  providing  a  sample  size  of  15.  The  length 
is  maintained  at  50  pm  while  the  total  number  of  voxels  used  to 
measure  the  sample  is  increased.  The  size  of  the  model  is  mea¬ 
sured  in  the  number  of  divisions  along  each  side  of  the  cube.  The 
physical  dimensions  of  the  particle  diameters  and  the  porosity’s 
characteristic  length  are  also  constant  for  each  size.  The  mean  (/) 
is  listed  and  the  standard  deviation  (S)  is  also  recorded.  The  same 
procedure  is  repeated  in  Table  3  for  the  CTE.  The  models  selected 
as  converged  are  shaded  in  light  grey. 

3.2.2.  RVE  size 

Testing  of  the  RVE  size  means  that  the  dimensions  of  the  sample 
size  will  change  while  the  physical  size  of  each  voxel  stays  constant. 
Table  4  shows  the  modulus  results,  and  the  RVE  size  is  varied  from 


Table  2 

Discretization  of  modulus  at  25  °C 


Model  # 

Number 

Division 

Length  (pm) 

d  (pm) 

X  (GPa) 

S 

%  Difference 

1 

15 

10 

50 

1 

84.97 

8.60 

- 

2 

15 

20 

50 

2 

78.33 

4.39 

7.82 

3 

15 

30 

50 

3 

73.78 

5.10 

5.81 

4 

15 

40 

50 

4 

71.25 

3.91 

3.43 

5 

15 

50 

50 

5 

69.29 

3.26 

2.76 

6 

15 

60 

50 

6 

69.51 

3.96 

0.32 
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Table  3 

Discretization  of  CTE  at  1000  °C 


Model  # 

Number 

Division 

Length  (pm) 

rf  41m) 

x  (°C  x  1 0‘) 

S 

%  Difference 

1 

15 

10 

50 

1 

14.04 

0.26 

2 

15 

20 

50 

2 

14.20 

0.14 

1.13 

3 

15 

30 

50 

3 

14.19 

0.09 

0.09 

Table  4 

RVE  Size  of  modulus  at  25  °C 


Model  # 

Number 

Division 

Length  (pm) 

(l  (|am) 

X  (GPa) 

S 

%  Difference 

7 

15 

10 

10 

5 

- 

- 

- 

8 

15 

20 

20 

5 

37.25 

9.12 

- 

9 

15 

30 

30 

5 

70.16 

5.83 

88.34 

10 

15 

40 

40 

5 

70.09 

4.33 

0.09 

11 

15 

50 

50 

5 

69.29 

3.26 

1.15 

12 

15 

60 

60 

5 

70.66 

2.45 

1.98 

Table  5 

RVE  Size  of  CTE  at  1000  °C 


Model  # 

Number 

Division 

Length  (pm) 

rf  (|im) 

X  ("C"1  x  10  <4 

5 

%Difference 

13 

15 

10 

16.7 

5 

- 

- 

- 

14 

15 

20 

33.3 

5 

14.10 

0.24 

- 

15 

15 

30 

50.0 

5 

14.19 

0.09 

0.63 

11 

9 

40 

66.7 

5 

14.23 

0.05 

0.27 

the  result  of  50  divisions  from  Table  2.  Table  5  repeats  the  procedure 
for  CTE. 

3.2.3.  Final  realizations 

Box  plots  of  the  mean  and  standard  deviation  of  the  modulus 
illustrate  the  differences  between  the  two  different  kinds  of  con¬ 
vergence.  In  Fig.  4a  the  mean  of  the  modulus  decreases  significantly 
as  the  voxel  length  is  decreased,  until  both  the  mean  and  standard 
deviation  do  not  vary  significantly.  However,  in  Fig.  4b  the  mean 
varies  little  after  the  first  sample,  while  the  standard  deviation 
continues  to  vary. 

The  measurements  of  CTE  showed  slightly  different  behavior  by 
being  far  less  dependent  on  voxel  size,  but  slightly  more  dependent 
on  RVE  size.  Yet,  since  thermal  expansion  is  primarily  a  constraint 
issue,  far  fewer  elements  are  needed  for  accurate  calculation  of  CTE 
coefficients. 

The  final  models  used  for  further  analysis  are  shown  in  Fig.  5a 
and  b.  A  quarter  section  of  the  model  is  removed  to  show  the 
microstructure  within  the  3D  reconstruction. 

3.3.  Computational  expense 

The  ability  to  select  the  minimum  model  needed  for  analysis  is 
significant  due  to  the  computational  expense.  Since  the  reconstruc¬ 
tions  are  three-dimensional  a  slight  increase  will  eventually  lead  to 
significant  increases  in  the  time  to  reconstruct  Ni-YSZ.  This  is  born 
out  in  Fig.  6  which  plots  the  clock  time  per  second  versus  the  model 
divisions.  However,  it  should  be  noted  that  computational  expense 
is  strictly  a  matter  of  total  elements  and  complexity  of  the  reference 
functions  used. 

3.4.  Material  properties 


3.4 A.  Modulus 

Based  on  the  discussion  above,  models  with  50  divisions  were 
selected  for  computing  the  elastic  modulus.  Fifty  different  realiza- 
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Nickel 


Nickel 

YSZ 


Pores 


YSZ 


Pores 


4.  Discussion 

4  A.  Realization  reconstruction 

The  simulated  annealing  method  efficiently  recreated  many 
different  realizations  of  the  given  Ni-YSZ  microstructure.  The  2- 
point  probability  function  also  proved  sufficient  in  recreating  the 
microstructure  so  long  as  each  of  the  three  phases  was  described. 
In  Fig.  1  the  plots  of  the  nickel  and  YSZ  match  almost  exactly 
the  analytical  functions  of  overlapping  spheres.  Additionally,  for 
all  realizations  listed  in  Tables  2-5  the  final  reconstructions  had 
probability  functions  identical  to  the  reference  functions. 

The  importance  of  using  the  pore  probability  function  is  clear  in 
Fig.  2.  The  function  matches  the  bounds  in  Eqs.  (2.6)  and  (2.7),  which 
means  the  pore  phase  behaves  in  a  random  manner.  The  short- 
range  order  in  the  pore  phase  is  also  a  good  fit  with  the  short-range 
order  in  the  combined  nickel  and  YSZ  phases.  This  is  reasonable; 
since  the  initial  realization  was  created  using  only  the  two  overlap¬ 
ping  sphere  functions,  the  remaining  pore  phase  would  have  been 
forced  in  a  controlled  distribution  around  the  spheres. 

The  lineal  chord  plot  in  Fig.  3  provides  limited  information  on 
the  clustering  behavior  of  the  microstructure.  The  majority  of  chord 
lengths  are  limited  to  the  diameter  size  for  all  three  phases.  The  YSZ 
phases  show  slightly  more  clustering  than  the  nickel  phase,  which  is 
reasonable  since  YSZ  has  a  higher  volume  fraction.  The  fact  that  the 
pore  phase,  which  has  the  highest  volume  fraction  in  the  material, 
shows  less  clustering  would  be  a  cause  of  the  short-range  order  in 
its  2-point  probability  function. 


Fig.  5eps.  Final  realizations  of  Ni-YSZ  for  the  (a)  CTE  model  with  30  divisions  and 
(b)  Young’s  modulus  model  with  50  divisions. 

tions  were  constructed  and  analyzed  in  each  coordinate  direction 
to  obtain  the  effective  elastic  modulus.  The  results  are  best  fit  to  a 
normal  distribution,  even  though  there  is  a  right  leaning  skew  as 
shown  in  Fig.  7. 

3.4.2.  Coefficient  of  thermal  expansion 

Since  the  CTE  had  very  small  standard  deviations,  the  FE  analysis 
focused  on  the  CTE  variation  due  to  change  in  temperature.  Fig.  8a 
shows  the  nickel,  Ni-YSZ,  and  YSZ  CTEs  for  0-1000  °C.  The  CTE  is 
linear  except  at  the  Curie  point  of  nickel  at  349  °C.  The  Curie  point 
behavior  is  zoomed  in  on  in  Fig.  8b. 


Fig.  7.  Statistical  variation  of  Young’s  modulus  for  50  total  realizations  measured 
in  three  directions  with  a  model  size  of  50  divisions  at  25  °C,  a  particle  diameter  of 
5  |jim,  and  a  volume  length  of  50  |xm. 


J.  Johnson,  J.  Qu  /  Journal  of  Power  Sources  181  (2008)  85-92 


91 


4.2.  Model  convergence 

The  two  different  convergence  tests  had  different  results 
depending  on  the  material  property  measured.  While  the  Young’s 
modulus  was  dependent  on  both  voxel  size  and  the  total  volume 
size,  the  CTE  was  insensitive  to  changes  in  either.  By  comparing  the 
discretization  and  RVE  results  plotted  in  Fig.  4,  different  conclusions 
can  be  drawn  about  the  FE  determination  of  the  Young’s  modulus. 
First,  as  voxel  size  decreases  the  modulus  will  be  more  accurately 
captured,  as  the  voxels  more  accurately  capture  the  microstruc¬ 
tures’  behavior.  Secondly,  once  a  sufficient  RVE  is  found,  a  continued 
increase  of  volume  will  affect  the  standard  deviation  of  the  modu¬ 
lus  while  the  mean  stays  constant.  This  makes  sense;  as  the  volume 
grows  it  behaves  in  an  increasingly  homogenous  manner.  Ther¬ 
mal  expansion,  as  a  constraint  problem,  was  very  insensitive  to 
discretization  error  and  was  much  more  dependant  on  RVE  size. 

4.3.  Material  properties 

Compared  to  the  Hashin-Shtrikman  upper  bounds  for  the 
Young’s  modulus,  these  results  are  very  reasonable  as  shown  in 
Fig.  9  [27].  Work  by  Radovic  et  al.  has  predicted  the  modulus 
results  for  40%  porosity  to  be  approximately  60  GPa,  which  is  22% 
lower  than  our  calculated  value  [26].  However,  since  the  real¬ 
ization  microstructure  is  generated  using  an  exact  description  of 
overlapping  spheres,  it  is  safe  to  assume  the  microstructures  are 
significantly  different. 

The  slight  right  leaning  skew  of  the  modulus  results  (Fig.  7)  sug¬ 
gests  that  the  modulus  tends  upward.  The  exact  reason  for  this  is 
unknown,  but  could  be  because  a  lower  limit  of  the  modulus  is 


Fig.  8.  CTE  of  Ni-YSZ  with  respect  to  temperature  for  (a)  0-1000  °C  and  (b)  around 
the  Curie  point. 


more  a  function  of  volume  fractions  than  microstructural  arrange¬ 
ment.  Our  lower  bound  of  60  GPa  does  compare  with  experimental 
results. 

The  FE  results  for  CTE  match  closely  experimental  results,  which 
range  from  10.5e-6  to  14e-6  for  100-950  °C  [28].  It  is  interesting  that 
CTE  measure  had  such  small  standard  deviations  for  every  model 
configuration  analyzed.  From  Tables  3  and  5  it  appears  that  CTE 
is  primarily  a  function  of  volume  fraction.  Also  the  CTE  requires  a 
sufficient  RVE  size  to  be  accurately  measured.  The  plots  of  CTE  for 
Ni-YSZ  (Fig.  8)  show  that  the  CTE  is  almost  equally  divided  between 
the  constituent  materials,  and  that  it  also  has  linear  behavior  except 
at  the  Curie  point. 

5.  Conclusions 

The  two  point  correlation  functions  for  overlapping  spheres 
were  used  to  digitally  recreate  the  material  nickel  zirconia  used 
for  SOFC  anodes.  Next  finite  element  analyses  were  used  to 
determine  the  Young’s  modulus  and  CTE  of  Ni-YSZ.  The  FE  mod¬ 
els  accurately  captured  the  microstructural  variation  within  the 
composite.  Study  of  multiple  realizations  showed  that  different 
microstructures  with  the  same  statistical  features  behave  in  a 
relatively  uniform  manner  after  convergence  is  reached.  To  that 
end,  the  reconstruction  method  can  be  used  to  create  models 
for  FE  analysis  to  study  the  stochastic  nature  of  metal-ceramic 
composites. 

Although  only  elastic  modulus  and  CTE  are  considered  in 
this  paper,  the  reconstruction  method  developed  here  can  also 
be  used  to  study  other  mechanical  and  thermal  behaviors  of 
multi-phase  composites,  such  as  fracture,  fatigue,  creep,  and 
effective  thermal  conductivity,  as  well  as  flow  properties  such 
as  permeability.  Results  of  these  studies  will  be  published 
elsewhere. 
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